#1.安装包
BiocManager::install("DESeq2")
# 加载包
library(DESeq2)   
library(ggplot2)   # 用于作图的包
#读入和处理数据
setwd("C:/Users/86183/Desktop/data for class3") #设置新路径
countData <- read.csv("mRNA_exprSet.csv")
rownames(countData)<-countData[,1]#蛋白质名设置为行名
countData<-countData[,-1]#删除前两列非数据列
 
#差异表达分析
condition=c(rep("0",123))#现将condition赋值为0
colData <- data.frame(row.names=colnames(countData), condition)#建表colData
for(i in 1:123)#建立循环，对condition分组
{
  if (substring(rownames(colData)[i],14,15)=="01")
    {colData$condition[i]="cancer"}
  else
   { colData$condition[i]="normal"}
}

dds <- DESeqDataSetFromMatrix(countData = countData, 
                              colData = colData, 
                              design = ~ condition)
dds2<-DESeq(dds)#差异表达分子
res<-results(dds2)#将结果用result()函数来获取
#查看结果
summary(res)
# res格式转化：用data.frame转化为表格形式
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
# 依次按照pvalue值log2FoldChange值进行排序
res1 <- res1[order(res1$pvalue, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
# 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
res1_up<- res1[which(res1$log2FoldChange >= 1 & res1$pvalue < 0.05),]      # 表达量显著上升的基因
res1_down<- res1[which(res1$log2FoldChange <= -1 & res1$pvalue < 0.05),]    # 表达量显著下降的基因
res1_total <- rbind(res1_up,res1_down)

#火山图可视化
genes<- res1
# 根据上调、下调、不变为基因添加颜色信息
genes$color <- ifelse(genes$padj<0.05 & abs(genes$log2FoldChange)>= 1,
                      ifelse(genes$log2FoldChange > 1,'up','down'),'normal')
color <- c(up="#F08080",down = "#00CED1",normal = "#969696")

jpeg("火山图.jpg", width = 800, height = 1000)
p <- ggplot(
  # 指定数据、映射、颜色
  genes, aes(log2FoldChange, -log10(padj), col = color)) +  
  geom_point(size=2) +
  theme_bw() +
  scale_color_manual(values = color) +
  # 辅助线
  labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) 
p
dev.off()

